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ABSTRACT: The development of tissue engineering hollow 
fiber bioreactors (HFB) requires the optimal design of the 
geometry and operation parameters of the system. This 
article provides a strategy for specifying operating condi- 
tions for the system based on mathematical models of 
oxygen delivery to the cell population. Analytical and 
numerical solutions of these models are developed based 
on Michaelis-Menten kinetics. Depending on the minimum 
oxygen concentration required to culture a functional cell 
population, together with the oxygen uptake kinetics, the 
strategy dictates the model needed to describe mass trans- 
port so that the operating conditions can be defined. If 
Cmin2>J*^m wc Capture oxygen uptake using zero-order 
kinetics and proceed analytically. This enables operating 
equations to be developed that allow the user to choose 
the medium flow rate, lumen length, and ECS depth to 
provide a prescribed value of c^i^. When Cmin ^^mi we use 
numerical techniques to solve full Michaelis-Menten 
kinetics and present operating data for the bioreactor. 
The strategy presented utilizes both analytical and numerical 
approaches and can be applied to any cell type with known 
oxygen transport properties and uptake kinetics. 
Biotechnol. Bioeng. 2011;108: 1450-1461. 
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Introduction 

HoUow fiber bioreactors (HFBs) are ideal for tissue 
engineering on a clinical scale because the large surface 
area to volume ratio wiU reduce the requirements of 
reagents, labor, and space: a hollow fiber system can be used 
to culture the same number of ceRs in 0.58 L as 1 m^ using 
standard flask culture techniques (EUis et al., 2005), and 
large cell numbers of up to 2 x 10** cell/mL can be obtained 
(Scragg, 1991). Knazek et al. (1972) were the first to report 
using a HFB for mammalian cell culture; since then the use 
of HFBs for mammalian cell expansion has become well 
documented (Tharakan and Chau, 1986) and several cell 
types have been cultured in HFBs including lymphocytes 
(Gramer and Poeschl, 2000; Gloeckner and Lemke, 2001), 
hepatocytes (Nyberg et al, 1994), and the osteogenic cell line 
560pZIPv.neo (Ellis and Chaudhuri, 2007). There is 
extensive understanding of HFB fluid dynamics and mass 
transport obtained from experimental and modeling 
studies, and a wealth of data on tissue physiology and cell 
metabolism in vivo and in vitro. For example, Abdullah et al. 
(2009) and Abdullah and Das (2007) have focused on high- 
density bone cell populations, whereas hepatocyte culture 
has provided a focus for bioartificial liver development 
through studies such as Hay et al. (2000), Kawazoe et al. 
(2006), Nyberg et al. (2005), Patzer (2004), Sielaff et al. 
(1997), Sullivan et al. (2007), and Wurm et al. (2009). 
Together these studies provide insight into the interaction 
between the cell environment and the fluid dynamics and 
mass transfer of nutrients across the membrane. Oxygen is 
recognized as the limiting nutrient with respect to growth of 
a cell population and has been the most widely modeled 
(although glucose has also been considered). The uptake of 
oxygen is usually modeled using Michaelis-Menten kinetics, 
which captures the dependence on the uptake rate on the 
underlying concentration. 
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As a consequence of the nonlinear nature of Michaelis- 
Menten kinetics, numerical solutions to the transport 
equations associated with HFBs are commonly seen in 
literature. These use full Michaelis-Menten; examples of 
finite difference methods include PUlarella and Zydney (1990), 
whereas examples of finite element methods include AbduUah 
and Das (2007), Chen and Palmer (2010), Das (2007), Sullivan 
et al. (2007, 2008), and Ye et al. (2006). Analytical approaches 
have also been used in literature for situations where 
Michaelis-Menten can be approximated by zero- or first- 
order kinetics. Example of zero-order kinetics are Piret and 
Cooney (1991), whereas examples of first-order kinetics are 
Jayaraman (1992) and Kim and Cooney (1976). Although 
Kim and Cooney ( 1976) use first-order kinetics, the ftmctional 
forms for the substrate concentrations that they determine are 
not dissimilar to those presented in this article. A good review 
of a range of transport models in HFBs is given by Brotherton 
and Chau (1996). 

To ensure the efficacy of HFB for clinical applications it is 
necessary to have information that allows accurate and 
correct operation of the HFB. This article presents a tool to 
select the modeling approach best suited to obtain cell type- 
specific operating data. As such, the approach presented 
here differs significantly from existing studies in the 
literature. First of all, previous studies have considered only 
analytical or numerical solutions in isolation. Here we use 
both approaches, and specify how to differentiate between the 
two based on cell data. Secondly, the analytical solutions that 
we present are based on zero-order kinetics and have not been 
reported previously in the literature. Finally, a strategy is 
outlined for providing operating data (specifically the lumen 
length, extra-capillary space (ECS) depth, and lumen flow 



rate) that ensure the oxygen concentration throughout a HFB 
is held above a prescribed tissue-specific minimum. When an 
analytical approach is applicable this data takes the form of 
operating equations that relate the underlying parameters; 
for the numerical approach operating data are presented 
graphically. This strategy enables a user to fix the geometry 
(e.g., lumen length, ECS depth) and operating conditions 
(e.g., lumen length) of the bioreactor to obtain their required 
cell culture environment. 



Theory 

Setup 

The fibers in a HFB fiber bundle are assumed to be Krogh 
cylinders, so that each fiber is identical and surrounded by 
an annulus of ECS containing a homogeneous distribution 
of cells (Krogh, 1918). The interstitial space between the 
Krogh cylinders is neglected as a modeling assumption. In 
this study, we consider transport in a single Krogh cylinder 
unit of a HFB bundle. This unit consists of a central lumen 
with a synthetic porous wall (referred to as the membrane), 
and surrounding ECS containing cells. Let z be the axial 
direction down the lumen centerline, starting at the lumen 
inlet (z = 0) with the lumen outlet denoted by z~L. We 
denote the radius of the lumen by d, the depth of the 
membrane by s and the depth of the ECS by I. Typical values 
are L=10cm, d=100|jLm, 5=20|xm, and Z = 600|jLm 
(Ye et al, 2007), although these should be varied as part 
of the bioreactor design process. A schematic of the setup is 
given in Figure 1. 




□ lun 



Membrane 



I I Extra-Capillary Space (ECS) 



Figure 1 . a schematic of the HFB setup. The left-hand schematic shows the structure of a fiber bundle, comprising seven Krogh cylinder units. The right-hand schematic 
shows a cross-section through an individual fiber, including the fluid velocity profile in the lumen. 
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Culture medium is pumped through the lumen at an 
imposed flowrate. There is no flow through the inlet to the 
membrane or ECS, so that fluid enters the system through 
the lumen only. Although this medium includes a mixture of 
solutes and proteins, we consider the transport of oxygen 
alone in this article. This is a widely adopted approach 
in the literature as oxygen is generally considered to be the 
rate-limiting nutrient, and reduces the complexity of the 
modeling process (Martin and Vermette, 2005; Piret and 
Cooney, 1991). Oxygen is transported by both advection 
(by the fluid) and diffusion in the lumen. Furthermore, 
oxygen diffuses through the membrane and ECS, where it is 
taken up by the cell population. In the analysis that follows 
we assume that the cell population is homogeneously 
distributed throughout the ECS, and neglect expansion of 
the cell population so that the parameters describing oxygen 
uptake are constant in time. 

Fluid flow in the lumen is described by Poiseuille's law 
whereas flow in the membrane and ECS is neglected (this is a 
common modeling assumption for small aspect ratio HFB 
when there is not a significant pressure drop across the 
membrane or ECS (Brotherton and Chau, 1996; Piret and 
Cooney, 1991)). We denote this fluid velocity in the lumen 
by u — 2U(l—r^/d^)e-, where U is the mean velocity 
(ms~^), r is the radial coordinate, and is the unit vector in 
the z-direction. The oxygen concentration and flux are 
denoted by c (molm~^) and J (molm~^s~'), respectively, 
with subscripts I, m, and e denoting the values in the lumen, 
membrane, and ECS, respectively. The oxygen fluxes are 

J; = c,u-AVq, J„, = -D,„Vc,„, ~D,Vce, (1) 

where D/, D„„ and Df, are the diffusion coefficients for 
oxygen in the lumen, wall, and ECS, respectively (all 
assumed constant, wdth units m^s~'). The lumen oxygen 
flux is comprised of advection due to the fluid velocity, 
together with diffusion; the membrane and ECS fluxes are 
comprised of diffusion only. The conservation equations 
for the concentration of oxygen in each of the regions are: 

— + V • J/ = 0 in the lumen, 
at 

+ V ■ J„, in the membrane, (2) 

of 

^ + "^ ■ie + R{ce) = 0 in the ECS, 
at 

where the reaction term RicJ captures the uptake of oxygen 
by the cells. We will assume Michaelis-Menten kinetics 
for this reaction term, so that 

It is necessary to prescribe boundary conditions on 
the internal and external boundaries of the bioreactor. On 



the lumen/membrane and membrane/ECS boundaries we 
prescribe continuity of concentration and flux, so that 

c; = c„, and J; • n = J„, • n 

(4) 

on the lumen/membrane boundary. 



Cm ~ Q and J„, • n = Je • n 

(5) 

on the membrane/ECS boundary, 

where n is the unit outward pointing normal to the relevant 
surface. Finally we prescribe the oxygen concentration as 
(molm~^) at the lumen inlet (where may be chosen to 
suit the application under consideration), and impose no 
flux of concentration out of the outer ECS boundary, 

Ci ~ Cin on z = 0, and Je • n = 0 

(6) 

on the outer ECS boundary. 

The assumption of no flux out of the outer boundary is 
analogous to a symmetry condition representation of a 
bundle of fibers. It compares directly to the Krogh cylinder 
approach used frequently in the literature. 

Next the solution of the model (2)-(6) is considered using 
numerical or analytical techniques. For both strategies a 
steady-state solution is sought and it is assumed that a 2D 
axisymmetric geometry is described by the radial coordinate 
r ~ ^Jx^^^Vy^ and the axial coordinate z. 

Analytical Approach 

To pursue an analytical approach, the system of equations 
given by (2)-(6) can be simplified with various assumptions. 
First of all the small aspect ratio of a fiber is exploited, 
defined by e = djL w 1 x 10"^ <C 1. It should be noted 
that whilst the lumen radius, d and fiber length, L can both 
be varied as part of the design process so that neither d nor L 
are fixed, s ^ 1 will be maintained throughout. 

It is not possible to make progress analytically using the 
nonlinear Michaelis-Menten reaction term given by (3). 
Therefore, we assume that q ^ K,n so that the reaction term 
_R(Cj,) can be approximated by V^^^^. This is an important 
assumption and means that predictions of the analytical 
model are only valid when the ECS oxygen concentration is 
much larger than the half-maximal oxygen concentration. 
As such, for cell types where the demand for oxygen is 
similar to, or smaller than, it will not be appropriate to 
use the analytical model (in this scenario a numerical 
approach should be used, as outlined later in the article). 

Finally the relative importance of advection and diffusion 
in the lumen is evaluated by considering the Peclet number, 
Pe= UL/Di. In fact it is the reduced Peclet number, 
Pe* = e^Pe = Ud^/LDi, that is critical for this system, as 
it also takes account of the small aspect ratio of the lumen (it 
is analogous to the reduced Reynolds number that was used to 
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characterize fluid transport for a similar study in Shipley et al. 
(2010)). A large reduced Peclet number indicates an 
advection-dominated regime, whereas a small reduced 
Peclet number indicate a difliision-dominated regime. 
Typically for this system C/wlcms~, LwlOcm, and 
DwlO~^m^s~\ giving Pe*wl so that advection and 
diffusion are both important in the lumen. It is assumed 
that Pe* = f^Pe is of order 1 in the analysis that follows. 
For the mathematical detail of the reduction of (2)-(6) based 
on the assumptions above, together with the solution of the 
resulting model, please refer to the Supplementary Material A. 

The outer radius of the lumen, membrane and ECS (each 
measured from the lumen centerline) are denoted by _R;, R„, 
and Re so that Ri^d, R^^d + s, and Rg = d + s + l. The 
following dimensionless parameters are also defined: 



and 



Pe* 



M 



d^V^ 



(7) 



which capture the key physical features of the system. As 
described above, Pe* is the reduced Peclet number and is 
assumed to be of order 1 . The parameter M represents the 
balance of oxygen consumption versus diffusion in the ECS, 
and can take a range of values depending on the relative 
importance of these effects. 

The analysis described above and in the Supplementary 
Material results in the following expressions for the oxygen 
concentration throughout the module: 



Q_ 



1 



d 



KummerMi 



V2 ■ 4 ' ' d^ 



G„exp 



2Pe*L kl 



iRi-Rl)ln- + Biz), 



(8) 



(9) 



M 



2R]\n 



+ B{z), 

where 

B(z) = 1 + y 



(10) 



+ X] £«exp KummerM + ^ ' ^ ' 



G„exp 



2Pe*L 



Fn 



(11) 



MDe 

2D,d^ 



(12) 



Here "KummerM (/x,v,x)" is the confluent hypergeo- 
metric function and is a solution of a specific differential 
equation, as described in the Supplementary Material A (and 
discussed in Abramowitz and Stegun, 1965). Further i„, £„, 
F„, and G„ for « = 0, . . . , oo are constants. The /,„ and E„ are 
the eigenvalues and normalization constants for the Sturm- 
Liouville problem associated with the system (2)-(4); these 
are constants independent of the geometry or cell 
population properties and are provided in the 
Supplementary Material B for « = 0, . . . , 49. By contrast 
F„ and G„ are coefficients in a Sturm-Liouville expansion of 
two different functions, and depend explicitly on the cell 
population properties (specifically the consumption rate of 
oxygen) and the geometry of the bioreactor (specifically the 
radius of the lumen and depths of the membrane and ECS). 

Although Equations (8)-(10) appear complex, the 
behavior that they describe is relatively straightforward to 
understand: the oxygen concentration in the lumen, 
membrane, and ECS depends on the radial distance from 
the lumen centerline. Each solution is also dependent on the 
distance down the lumen centerline, z, as a consequence of 
advection in the lumen. This is transmitted into the 
membrane and ECS regions through the function B(z), 
which is the lumen concentration value on the lumen wall 
(i.e., the solution in (8) when r^Ri = d). This function B(z) 
reveals that the concentration decays exponentially down 
the lumen fi"om a maximum value at the inlet z = 0. The 
remaining terms in the solution for c,„ and in (9) and (10) 
describe the radial decay of the oxygen concentration from 
the outer surface of the membrane as a consequence of 
oxygen uptake by the cells in the ECS. 

Through cell-specific design criteria, we must design the 
bioreactor to ensure that the oxygen concentration exceeds a 
prescribed minimum throughout the bioreactor. This 
minimum oxygen concentration will be achieved at the 
furthest distance from the inlet, that is, when r = Ri, and 
z = L. Denoting this minimum value by c^iny the analytical 
method gives the following expression for c^in, in terms of 
experimentally controlled and cell-specific parameters: 



M 
AcP 



Ri-Ri 



2D, 

Drii 



iRi-R:)^^ 



R„ 



-2Rlln 



R.. 



B(L). 



(13) 



Numerical Approach 

For the analytical approach, the full system given by (2)-(6) 
is solved using finite element method package "COMSOL 
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Multiphysics 3.5a"^ to evaluate the dependence of the 
oxygen concentration on the underlying parameters. The 
numerical approach is valid for all concentration values; 
however, the full system of equations must be solved 
iteratively each time. This is a computationally intensive 
process and does not provide operating equations that 
describe the dependence of the minimum oxygen concen- 
tration on the underlying parameters. Therefore, the 
numerical approach will be used when the analytical 
approach is not valid, that is, when c^nm/^K,,,. The mesh 
used for the results in this article consists of approximately 
7,000 finite elements (and refining the mesh to 29,312 
elements did not change the results to three significant 
figures). 

Results and Discussion 

The analytical and numerical methodologies outlined in 
the Theory Section will be used to outline a strategy for 
developing cell-specific operating criteria for the bioreactor. 
These criteria will then be tested for specific cell types. 



Strategy for Developing Optimal Operating Conditions 

To develop operating conditions, it is necessary to 
understand how the minimum oxygen concentration 
depends on the geometrical properties of the bioreactor, 
together with the parameters that can be controlled 
experimentally. Through this understanding, the HFB 
can be designed to optimally grow cells of a particular 
type. 

Once a cell type and seeding density are chosen the 
following parameters are fixed: 

(1) The maximal oxygen consumption rate, V^^,^. 

(2) The half-maximal oxygen concentration, K^. 

(3) The diffusivity of oxygen in the ECS, D^. 

The diffusivity of oxygen in the lumen and membrane (D; 
and D„„ respectively) are known from the literature or 
experiments. The outer radii of the lumen and membrane 
{Ri and R^, respectively) are fixed, and there are specific 
values of the lumen inlet concentration c-^ and minimum 
oxygen concentration c,nin that must be achieved. So, the 
bioreactor design parameters that are left to be determined 
are: 

(1) The depth of the ECS, I (which determines R^). 

(2) The length of the lumen, L. 



'Developed and distributed by COMSOL, Inc. Full details available online at 
http://www.comsol.com/. 



Finally, the mean inlet flow rate U can be controlled by 
fixing the volumetric flow rate on the pump used to deliver 
fluid to the bioreactor. 

If c^in^ K„, the analytical approach is valid and the 
results from the Analytical Approach Section can be used to 
fix 1, L, and U; however, if Cmin ^ K„, the analytical approach 
is not valid, and the numerical method must instead be used. 
These two approaches are detailed below. 

Ciain^K,!! The operating conditions are specified for 
the bioreactor using Equation (13) for the minimum 
oxygen concentration. When the parameters described 
above are fixed, only .R^ and the ratio U/L (through Pe*) 
can be determined independently using the analytical 
approach. Two cases will be considered: 

(1) The outer radius of the ECS, R^ is fixed, and so Pe* can 
be determined. 

(2) The ratio U/L (and therefore Pe*) is fixed, and so R,. can 
be determined. 



For the first case it is assumed that the outer radius of the 
ECS is fixed so that R^ (and thus y, F„, and G„ for 
n = 0, . . . , oo) is known. In this case, (13) can be written as 
the following operating equation for c^i^ in terms of the 
reduced Peclet number Pe*: 

^ = A + £B„expf--^)+C„, (14) 

where 

+ 1 + 7, 

(15) 

B„ = _E„G„exp KummerM Q + ^ ' ^ ' ' 

C„ — exp I — I KummerM | - H — - ,\, —X,, ] , 

(16) 

are all fixed constants. Given the values of these constants, 
(14) can be used to determine the value of Pe* (and hence 
the ratio U/L) that provides the required value of c^i^ (note 
that it is this ratio rather that the individual values of U and 
L that influence the minimum oxygen concentration). 
Equation (14) shows that the minimum oxygen value c^i^ 
decreases exponentially as the lumen length L increases, or 
the lumen velocity U decreases. This means that for a lower 
Cmin requirement, a smaller flow velocity and longer fiber can 
be used. 
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Next it is assumed that the ratio UIL is prescribed so that 
Pe* is given. Now the operating equation for Cmin in terms of 
the ECS depth is: 



Cell Types and Parameter Values 

The parameters that will be kept fixed in our investigation 



where 



K 



M 
AcT- 



D,„ \ a 



1 



M [ fR„ 



1- 



De \ 2D, 



D„ 



D, 



(17) 



(18) 



H„ = E„exp KummerM ^ + ^ > 1 > -^«^ exp (- ^ 

En A'A A ^« A 
Jn — ^ exp — KummerM - H , 1 , —X„ , 



(19) 



(1) The oxygen diffusivities D;=3xlO 'm^s \ 
D„ = 3 X 10"'°m^s"S and D^ = 6 x 10"''m^s"^ (Ye 
et al., 2006). 

(2) The lumen radius _R;= ii= 100 |jLm and the depth of 
the membrane s = 20|xm (so that _Rm=120p-m) (Ye 
et al, 2006). 

(3) The inlet oxygen concentration will be fixed for each 
individual cell type. 



The kinetic data (i.e., Vj^a^ and i<C,„) for a range of cell 
types, sourced from combined modeling and experimental 
studies in the literature, are shown in Table I. For 
cardiomyocytes, hepatocytes, and pancreatic cells we fix 
Cin = 0.22molm~ (as is standard for culture medium 
Piret and Cooney, 1991). However, chondrogenic differ- 
entiation is limited when the oxygen concentration exceeds 
approximately O.lmolm"^ (Lund-Olesen, 1970; Treuhaft 
and McCarty, 1971); therefore c^^ = 0.1 molm~^ is used for 
chondrocytes. 



are all fixed constants. Given the values of these constants, 
(17) depends on R,. through the explicit appearance of R^. in 
(17) as well as G„ and F„ for « = 0, . . . , oo. For a given value 
of Cjnin, (17) can therefore be solved numerically to 
determine R,.. 

Cmin In this scenario the numerical approach 

will be used, as outlined in Numerical Approach Section. 



Validation of Analytical and Numerical Approaches 

The analytical approach is a reduction of the full model 
given by (2)-(6) and therefore should be validated. This 
validation could be performed against experimental data; 
however, this data is difficult to coUect accurately and is not 



Table I. Oxygen uptake and culturing data for a range of cell types. 





^max 


Km 


Cell density 
(cells m"') 


^min 






Cell type 


(mol m^' s^') 


(mol m"') 


(molm"^) 


(molm"') 


Source 


Neonatal rat cardiomyocytes 


2.64 X 10"^' 


6.9 X 10"^' 


10'^ 


8 X 10"^ 
6 X 10"' 


0.22 


Radisic et al. (2005) 
Carrier et al. (1999) 


Primary rat hepatocytes 


1.76 X 10"'' 


6.24 X 10"^' 


1.25 X lO'' 


2.1 X 10"^ 


0.22 


Sullivan et al. (2007) 
Consolo et al. (2008) 


Pancreatic pTC3 cells 


6.37 X 10"^' 


1.0 X 10"^ 


2.8 X lO" 


1.46 X 10"^ 


0.22 


Tziampazis and Sambanis (1995) 
Stabler et al. (2009) 


Bovine chondrocytes 


4.8 X 10"^ 


5.0 X 10"' 


1.4 X 10" 


1 X 10"^ 
1.32 X 10"^ 
2.2 X 10"' 


0.1 


IVIalda et al. (2004) 
Obradovic et al. (1999, 2000) 
Fermor et al. (2007) 



For a description of the various minimum oxygen concentrations, please refer to the main text. The V„„ value for neonatal rat cardiomyocytes and 
primary rat hepatocytes have been multiplied by a cell volume fraction of 0.3, as per the modeling in Sullivan et al. (2007). For the pancreatic cells it has also 
been assumed that each cell has a 10 |j.m diameter. 

Note: For neonatal rat cardiomyocytes two values are listed. It has been observed that cardiac constructs cultivated in perfusion at oxygen concentrations of 
~80 (jlM exhibit weaker presence of cardiac markers and poorer organization of contractile apparatus compared with constructs cultivated at oxygen 
concentrations of ~200 (jlM Carrier et al. (1999); this explains the first value. The second value (6 p-M) is a typical hypoxia value (Radisic et al., 2005). The 
value for primary rat hepatocytes is based the critical threshold value of 10 mmHg quoted in the literature Consolo et al. (2008) (and transferred from a partial 
pressure into a concentration using Henry's law with an oxygen solubility value of 2.08mmolm"' mmHg). For pancreatic pTC3 cells, published experiments 
found that oxygen tensions above 7 mmHg were required for the cells to retain their secretory capacity Stabler et al. (2009); using Henry's law gives the value in 
Table I. Finally, a range of minimum oxygen concentrations are presented for articular cartilage in the literature. In Obradovic et al. (1999), it is hypothesized 
that articular cartilage is exposed to a minimum oxygen concentration in the range 0.01 mol m"' to 0.08 mol m"' in vivo, where lower oxygen concentrations 
are not detrimental to chondrocyte viability but can impact synthesis of extracellular matrix; this explains the first c„,i„ value in Table I. In Fermor et al. (2007), 
it is reported that the superficial zone of articular cartilage exists at above approximately 6% oxygen concentration, whereas the deep zone exists at <1%; this 
explains the final two c„i„ values of Table I. 
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- - z=L: Numerics 

— z=L: Analytics 

Oil!^ . 1 1 i 1 

0 50 100 150 200 250 

Radial distance r (u m) 

Figure 2. comparison of the analytical and numerical approaches. The graph 
shows the radial oxygen concentration profiles for primary rat hepatocytes (see 
Table I) at fixed values of z, using both the analytical and numerical techniques. The 
fixed parameters are 0= 1 x 10"^ ms"', /. = 10 cm, and /?;, = 220 jim. 



presented in sufficient detail in the hterature. Given that the 
numerical approach is valid for all concentration values and 
solves the full model (2)-(6), it is appropriate to validate 
results of the analytical model against numerical solutions. 
This comparison is shown in Figure 2, where radial oxygen 
concentration profiles are shown (at fixed values z = 0, L/3, 
2L/3, L) for the primary rat hepatocyte data in Table I 
(Sullivan et al., 2007) when U= 1 x 10"^m/s, L= 10 cm, 
Re = 220 ixm. For the analytical solutions, all sums have been 
truncated at 50 terms, that is, n = 49, for ease of 
computation. The agreement between the analytical and 
numerical results is very strong, although it becomes weaker 
as the concentrations decrease. The lowest concentration 
value is at the ECS outlet (when r^R^ and z^L); here both 
the analytical and numerical concentration values are 
0.12molm~^ to two decimal places, with a percentage 
difference of 3.39% (which is within experimental error). 

Analytical and Numerical Results 

It must first be decided whether to use the analytical or 
numerical strategy to provide operating data. Table II 
provides a summary of this decision making process. Data 
on Cjnin and are provided for each cell type, together 
with the value of the ratio c^iJK,,,. The analytical model 
is vahd when c^in 3> K„j; here we choose a value of the 
ratio c^iJK^ = 2 as the critical value so that if c^iJK^ > 2 
the analytical model is used, whereas if c^iJK^ > 2 the 
numerical approach is used. Different critical values of c^^J 
could certainly be implemented, even on a cell-specific 
basis. The errors associated with using c^iJK„, = 2 as the 
critical value are within the bounds of experimental error, 
and the errors associated with other modeling assumptions 
(e.g., the Krogh cylinder approximation). On this basis, the 
analytical model is appropriate for the cardiomyocytes 
(cmm=8x 10~^molm~^), hepatocytes, and chondrocytes 
(Cmm= 1-32 X 10~^molm~^^), whereas the numerical model 
is used for the remaining examples in Table II. We present 



Table II. Use of the analytical or numerical models. If c^^JK„, > 2, the 
analytical model is used; otherwise the numerical model is used. 



Analytical Numerical 
Cell type c^ijci^ KiJc^a c^iJK^ model model 



Neonatal rat 


0.36 


0.031 


11.6 




X 


cardiomyocytes 














0.027 


0.031 


0.87 


X 




Primary rat 


0.095 


0.028 


3.4 


1^ 


X 


hepatocytes 












Pancreatic 


0.066 


0.045 


1.5 


X 


\^ 


PTC3 ceUs 












Bovine 


0.1 


0.05 


2.0 


X 




chondrocytes 














0.13 


0.05 


2.64 


1^ 


X 




2.2 X 10"^ 


0.05 


0.44 


X 




data for the 


extreme 


cases 


of high 


and low 


oxygen 



requirements, that is, cardiomyocytes and chondrocytes, 
respectively. 

Figures 3a-c and 5a-c show the variation in c^yjc-^^ (with 
fixed) as a function of l/_Pe* for fixed R^, as described by 
operating equation (14), for the cardiomyocytes and 
chondrocytes, respectively. As would be anticipated, c^^^ 
is largest for low values of l/Pe*, corresponding to either a 
large lumen velocity [/or shorter lumen length L (a larger U 
ensures increased delivery of oxygen to the cells through 
advection, whereas a shorter lumen length decreases the 
distance of the furthermost cells from the oxygen source). 
For each value of R^ the maximum variation in is of 

size 10~ , indicating that c^i„ is only weakly sensitive to the 
value of Pe*. For Pe* < 2 (i.e., l/Pe* > 0.5) c„,in is virtually 
constant, indicating a linear relationship between the values 
of 17 and L required to achieve a chosen value of c^nin. 

Given that this linear relationship is representative of the 
low Cjnin regime, it is mimicked by the numerical results of 
Figures 4 and 6. These figures show how the critical lumen 
length, Lcrit say, required to satisfy the minimum oxygen 
concentrations of Table I varies as a function of the lumen 
velocity U. For each cell type four different values of 
Re were tested, each of which demonstrates a linear 
relationship between Lerit and U (with correlation factor 
0.99). These figures can be used to read off a required 
Lcrit and U value to satisfy the minimum oxygen 
requirements summarized in Table I. For example, for 
the cardiomyocytes with i?^ = 270 |xm with a lumen length 
of 10 cm, a lumen flow velocity of [/w9 x 10~''ms~' wiU 
ensure c>6 x 10~ molm~^ throughout the module. 
In contrast. Figures 3d and 5d show the variation in c^iJ 
(with fixed in each case) as a function of R^ for fixed 
Pe*, as described by operating equation (17), for the 
cardiomyocytes and chondrocytes, respectively. As antici- 
pated, Cmi„ decreases as the ECS depth (i.e., R^.) increases. 
The rate of this decay is heavily dependent on the uptake rate 
of oxygen by the cell population, V„^^„ which is largest for 
the cardiomyocytes, and lowest for the chondrocytes. For 
example, a value of Cmin/cin = 0-2 is sustained by the 
cardiomyocytes, hepatocytes, and chondrocytes when 
R^ w 267 and 720 |jLm, respectively. 
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Case Study: Cells With a High Oxygen Requirement 

Consider a HFB for culturing cardiomyocytes with the 
following known parameters: Crni„ = 8xl0~ molm~ , 
Cin = 0.22molm~ , and c^^iJK„,= 11.6, together with the 
ECS depth fixed at 95 |jLm so that i^p = 215|xm. For HFB 
operation it is necessary to specify the inlet flowrate and 




"a 4 6 8 10 12 14 

Mean Lumen Velocity U (m s~^) xlO^ 

Fiyure 4. Numerical results for the neonatal rat cardiomyocytes that show the 
relationship between and t/when c„\„ = 6 x lO^^mol m"^ and c^JK^= 0.87 are 
held fixed (arrow in direction of fl, decreasing). [Color figure can be seen in the online 
version of this article, available at http://wileyonlinelibrary.com/bitl 



fiber length to maintain the oxygen concentration above this 
minimum. Since c^^^^JK^ > 2 and the ECS depth is fixed, the 
analytical approach (operating equation 14) will be used to 
determine the value of Pe* (and corresponding possible 
values of L and U) that achieves = 0.36. For this 

scenario, A = 0.80 and the values of 5„ and C„ for 
n = 0, ... ,49 are given in the Supplementary Material C. 
Solving operating equation (14) yields Pe* = 0.2 so that the 
ratio U/L = 6 x 10~^. Any values of U and L that satisfy this 
ratio will ensure c> 8 x 10~^molm~^^ throughout the 
construct; two example values are L = 0. 1 m and 
(7=6 X 10"^ms"\ 

By comparison, suppose the lumen flow velocity is fixed 
at U= 1 X 10~^ms~' and the lumen length at L = 0.1 m so 
that Pe* = 1/3. Then operating equation (17) can be used to 
determine the ECS depth that achieves c^iJcin — 0.36. For 
this scenario, i<:= 1.01 and Q= 1.18 x 10*m"^ and the 
values of H„ and /„ for n = 0, . . . , 49 are given in the 
Supplementary Material C. Solving operating equation (17) 
with Cjnin/cin = 0.36 uow givcs — 212.8 |jLm so that the ECS 
depth is 92.8 |xm. 

Case Study: Cells With a Low Oxygen Requirement 

Next consider a HFB for culturing chondrocytes with the 
following known parameters: Crnin = 2.2 x 10~ molm~ , 
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Cin = 0.1 molm~^, and Cmin/^T,,, = 0.44, together with the 
ECS depth fixed at 150|jLm so that _R^ = 270|xm. It is 
necessary to specify the inlet flowrate and fiber length to 
maintain the oxygen concentration above this minimum. 
Since c^-^JK,,, < 2 and the ECS depth is fixed, the numerical 
approach wiU be used to determine possible values of L and 
U that achieve c^^Jcin^ 2.2 x 10~^. For this scenario, we 
refer to Figure 6. Any values of L and U that lie on the blue 
line (_Re= 270 |JLm) are appropriate: an example is 
17= 2.8 X 10~''ms~^ and L = 0.2m.By comparison, sup- 
pose the lumen flow velocity is fixed at U= 3 x 10~''ms~^ 
and lumen length L = 0.5. The red line of Figure 6 dictates 
that Re — 170 |xm should be imposed in this case. 



Discussion 

The strategy that has been outlined enables mathematical 
modeling techniques to inform bioreactor design based on 
the oxygen requirements of the cell type. Two different 
modeling approaches were employed to provide design 
and operating data that ensure the oxygen concentration 
throughout a HFB is held above a prescribed tissue-specific 
minimum value, Cmin that ensures the grovrth of a functional 
cell population. When c^^^ 3> -fT,,, (the half-maximal oxygen 
concentration), oxygen uptake by the cell population was 



captured using zero-order kinetics, and operating equations 
were derived analytically. These operating equations provide 
insight into the relationship between the minimum 
oxygen concentration and the geometrical properties of 
the bioreactor, together with the operational parameters 
(such as inlet oxygen concentration and flow rate) than can 
be controlled by the user. A case study was presented that 
demonstrated how to use these operating equations for cell 
types with a high oxygen requirement. However, an 
analytical approach is not valid when c^nm/^K,,,. In this 
case, full Michaelis-Menten kinetics must be solved in the 
ECS using a numerical approach. This was achieved using 
the finite elements package "COMSOL Multiphysics," and 
operating data on the relationship between lumen length 
and flow rate required to achieve a specific minimum 
oxygen concentration value were presented. This approach 
has the advantage of being valid for aU concentration values; 
however, the full system of equations must be solved 
iteratively each time and this is a computationally intensive 
process. A case study was presented that demonstrated 
how to use these operating equations for cell types with a low 
oxygen requirement. 

Previous studies into the modeling of tissue engineering 
bioreactors have focused on either numerical or analytical 
approaches (under various simplifying assumptions) in 
isolation. For example, AbduUah and Das (2007), Chen 
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Figure 6. Numerical results for the bovine chondrocytes that show the relationship between /.critand Ufortwo different minimum oxygen requirements (arrows in direction of 
flg decreasing), a: [Cmi„= 1 x lO^^mol m"^ and c„i„IK„, = 2.0 held fixed], (b) [c„i„ = 2.2 x lO^^mol m"^ and c„ JK„= 0.44 held fixed]. [Color figure can be seen in the online version 
of this article, available at http://wileyonlinelibrary.com/bit] 



and Palmer (2010), Das (2007), Pillarella and Zydney 
(1990), Sullivan et al. (2007, 2008), and Ye et al. (2006) 
employed various numerical techniques to solve full 
Michaelis-Menten kinetics for individual cell types in 
HFBs. By comparison, analytical approaches such as Piret 
and Cooney (1991), Jayaraman (1992), and Kim and 
Cooney (1976) have been used to approximate Michaelis- 



Menten by zero- or first-order kinetics. However, zero-order 
kinetics have not previously been used to determine 
operating equations, whilst first-order kinetics are only 
valid when the substrate concentration is smaller than the 
half-maximal substrate concentration, K^- This is not 
appropriate in the development of oxygen-based operating 
equations for the use of HFB for tissue engineering, where 
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the oxygen concentration must typically be maintained 
above K„ to ensure the growth of functional tissue. While 
these are all valid and workable models, they have not 
previously been integrated to provide a strategy that could 
be applied to any cell type to stipulate bioreactor design and 
operation. 

Conclusion 

A strategy has been developed for modeling oxygen kinetics 
in tissue engineering HFB. The strategy allows operating 
parameters to be specified that ensure the oxygen 
concentration is maintained above a prescribed minimum 
throughout the HFB. The strategy dictates that the 
appropriate approach is based on whether the Michaelis- 
Menten kinetics can be reduced to zero-order; in the case 
of high oxygen requirements zero-order kinetics is appro- 
priate and so the analytical approach is used. In the case 
of low oxygen requirements it is necessary to use fuU 
Michaelis-Menten kinetics and so a numerical approach is 
required. As such, the strategy developed here can be used 
for any cell type to specify operating parameters. 



M dimensionless constant that represents the balance 

of oxygen consumption versus diffusion in the 
ECS 

-y algebraically convenient parameter that depends 

on Rg 

A„ eigenvalues of Sturm-Liouville problem for 

n = 0, . . . , oo 

E„ normalization constants of Sturm-Liouville pro- 

blem for fi = 0, . . . , oo 

F„, G„ Sturm-Liouville expansion constants for 

« = 0, . . . , oo 

B{z) dimensionless oxygen concentration on the 

lumen wall 

A, B„, C„, K, Q, H„, J„ constants associated with the analytic operating 
equations 
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Nomenclature 



d radius of the lumen (m) 

s depth of the lumen wall (m) 

/ depth of the ECS (m) 

L length of a single module (m) 

L^t critical length required to satisfy a minimum 

oxygen requirement (m) 

z axial length coordinate down the lumen 

r radial coordinate 

u fluid velocity vector (ms^') 

U Mean velocity in the lumen (ms^') 

unit vector in the z-direction 

c oxygen concentration (molm^'') 

J oxygen flux (molm^^ s^') 

U velocity scale (ms^') 

Di oxygen diffusion coefficient in the lumen (m^ s^ ') 

D„ oxygen diffusion coefficient in the wall (m^s^') 

Dg oxygen diffusion coefficient in the ECS (m^ s^') 

R uptake rate of oxygen (molm^^s^') 

V„„ Maximal oxygen consumption rate (mol m^' s^ ') 

K„ half-maximal oxygen concentration (mol m^^) 

n unit outward pointing normal to a surface 

c-^ fixed oxygen concentration at the lumen inlet 
(mol m^^) 

dnin minimum oxygen concentration in the HFB 
(molm^^) 

e aspect ratio of the lumen 

Pe axial Peclet number in the lumen 

Pe* reduced Peclet number in the lumen 

Ri outer lumen radius (m) 

R„ outer membrane radius (m) 

Rg outer ECS radius (m) 
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